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ABSTRACT 

Conservative numerical schemes for general relativistic magnetohydrodynamics 
(GRMHD) require a method for transforming between “conserved” variables such as 
momentum and energy density and “primitive” variables such as rest-mass density, in¬ 
ternal energy, and components of the four-velocity. The forward transformation (primi¬ 
tive to conserved) has a closed-form solution, but the inverse transformation (conserved 
to primitive) requires the solution of a set of five nonlinear equations. Here we discuss 
the mathematical properties of the inverse transformation and present six numerical 
methods for performing the inversion. The first method solves the full set of five non¬ 
linear equations directly using a Newton-Raphson scheme and a guess from the previous 
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timestep. The other methods reduce the five nonlinear equations to either one or two 
nonlinear equations that are solved numerically. Comparisons between the methods 
are made using a survey over phase space, a two-dimensional explosion problem, and 
a general relativistic MHD accretion disk simulation. The run-time of the methods is 
also examined. Code implementing the schemes is available for download on the web. 

Subject headings: hydrodynamics — methods; numerical — MHD 


1. Introduction 

It is commonly thought that many astrophysical systems contain relativistic plasmas with 
a dynamically significant magnetic field. Examples include accreting black holes in black hole 
binaries, galactic nuclei, gamma-ray bursts, the cores of massive stars undergoing core collapse, 
isolated neutron stars, and neutron stars in binary systems. 

As a result, there is currently considerable interest in numerical methods for integrating the 
equations of general relativistic magnetohydrodynamics (GRMHD). Within the last few years about 
six GRMHD schemes have been deployed (Komissarov 2005; Koide, Shibata, & Kudoh 1999; Gam- 
mie, McKinney, & Toth 2003; De Villiers & Hawley 2003; Fragile 2005; Duez et al. 2005; Anton et 
al. 2005). 

Some of these authors (Komissarov 2005; Koide, Shibata, & Kudoh 1999; Gammie, McKinney, 
& Toth 2003; Duez et al. 2005; Anton et al. 2005) have adopted a conservative scheme. This means 
that the integrated equations are of the form 

9tU(P) = -d,F*(P) + S(P) (1) 

Here U is a vector of “conserved” variables, such as particle number density, or energy or momentum 
density in the coordinate frame, the F* are the fluxes, and S is a vector of source terms that do 
not involve derivatives of P and therefore do not affect the characteristic structure of the system. 
U is conserved in the sense that, if S = 0, the rate of change of the integral of U over the volume 
depends only on fluxes at the boundaries, by the divergence theorem. The vector P is composed 
of “primitive” variables such as rest-mass density, internal energy density, velocity components, 
and magnetic field components. The fluxes and conserved quantities depend on P. Gonservative 
numerical schemes advance U, then depending on the order of the scheme, calculate P(U) once or 
twice per timestep. 

In nonrelativistic conservative MHD schemes the conserved quantities are trivially related to 
the primitive variables; both the forward transformation P ^ U and the inverse transformation 
U ^ P have a closed-form solution. In GRMHD (or even SRMHD) U(P) is a complicated, 
nonlinear relation. The inverse transformation has no closed-form solution and must be performed 
numerically. This numerical operation must be robust, accurate, and fast—it is at the heart of all 
conservative GRMHD schemes. 
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In this paper we investigate several schemes for the inversion P(U) and test each in an ax- 
isymmetric simulation of accretion onto a rotating black hole. § 2 covers definitions and notational 
matters. § 3 describes five distinct formulations of the algebraic equations to be solved numerically, 
and § 4 describes how each method is implemented numerically. § 5 describes the performance of 
these methods in the context of a survey over a range of primitive variable values and in two typical 
applications. § 6 summarizes our findings and contains a guide to our results for those wishing to 
make their own implementation. We will assume throughout that G = c = 1. 


2. Definitions 

Throughout this paper we follow standard notation (Misner, Thorne, & Wheeler 1970). We 
work in a coordinate basis with metric components g^i, and independent variables t, xi,X 2 ,X 3 . The 
normal observer’s four-velocity is = (— 0 , 0 ,0,0) in this coordinate basis, where c? = is 

the square of the lapse. 

The fluid is described by its four-velocity u^, rest-mass density po, internal energy per unit 
(proper) volume u, and pressure p. The electromagnetic field is described by the field tensor , 
which is antisymmetric, and its dual 

( 2 ) 

The field tensor has six degrees of freedom (three components each for E and B), but three are 
eliminated by the ideal MHD condition 


= 0 (3) 

which states that the Lorentz force vanishes in the rest-frame of the fluid. It is convenient to 
describe the field using the magnetic field four-vector 

= -nJ’F^''. (4) 


Notice that B^n^ = 0 and the vector B differs from the magnetic field variables 


= *F 


it 


a 


( 5 ) 


that we have used in earlier papers in this series (Gammie, McKinney, & Toth 2003; Gammie, 
Shapiro, & McKinney 2004; McKinney & Gammie 2004). It is also useful to define the projection 
tensors 

hfMu ~ 9^iu T (6) 

which projects into a space normal to the fluid four-velocity and 


jijtv — 9^iv T 


( 7 ) 
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which projects into the space normal to the normal observer n^. 

The governing equations for GRMHD are conservation of stress-energy, 

= 0 ( 8 ) 

conservation of particle number, 

= 0 , 

and the homogeneous Maxwell equations 

= 0 . 

Often, we will assume a T-law gas, with equation of state 

p = (T — 1) u. (11) 

Some (but not all) of the schemes described here do not work for other equations of state. We will 
note where our assumed equation of state is essential. 

The GRMHD stress-energy tensor is 

= (u; + b^) + (^ + y) 9^’' " ' (1^) 

Here 

w = po+p + u, (13) 

^2 = and 

(14) 

7 

where 

7 = -n^u^ (15) 

is the Lorentz factor of the flow as measured in the normal observer frame; in our coordinate basis 
7 = au*. 

There are many possible choices of conserved and primitive variables. For definiteness, we will 
make a specific choice that is nearly identical to that used in the HARM code (Gammie, McKinney, 
& Toth 2003). Once an inversion scheme is found for this particular set many other choices can be 
obtained by simple algebraic transformations. 

For primitive variables, we use po,u,B^, and u* = j'^ (u* = 0). These velocities are the 
projection of the plasma four-velocity into the space perpendicular to n^. They can be rewritten 
ip = or, in the language of the ADM (Arnowitt-Deser-Misner) formalism -u* = u^+'yf3^ ja, 

where /3* = is the “shift”. The iPs are numerically convenient because they range from —oo 
to oo. 


(9) 

( 10 ) 
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For conserved variables, it is helpful to first write out the basic equations in full. The equations 
of energy-momentum conservation are 

dt + d, ( 16 ) 

where g = Det( 5 '^i,). The associated conserved variables are It is convenient to convert 

these to 

= -n,T% = aT\, ( 17 ) 

which is the energy-momentum density in the normal observer frame. 

The Maxwell equations yield three evolutionary equations 

dt = -dj [^/^ [Vu" - Bw>)] ( 18 ) 

and the constraint 

dt i^gB^) = 0 , ( 19 ) 

which does not concern us here. 

The particle number conservation equation is 

dt [yf^PoU^) = -dj {yf^PoU^) ■ ( 20 ) 

The conserved variable for this equation is y/—gpoU^. It is convenient to convert this to 

D = -poTi^u^ = Poau^ = 7/3o • (21) 

This is the density measured in the normal observer frame. 

To sum up, the eight conserved variables are Q^, D, and We are now in a position to give 
explicit expressions for the forward transformation U(P). 

First, we calculate p = {T — 1) u, w = po + u + p, j = yjl + gijii^id, = ( 7 /a, — ajg^^), 

and Then 

D = iPo, ( 22 ) 

QtJi = l{w + b^) Ufi- {p + &^/ 2 ) n/, + {uyb'') b^, ( 23 ) 

and of course the B'^s are both primitive and conserved variables. Along the way, it is computa¬ 
tionally efficient to use the identities 

62 = J- \b^ + {B^^u^f] (24) 

T L J 

and 

n^6^ = -u^B^. (25) 

The next section will show how these equations are used to perform the inverse transformation. 
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3. Inversion Schemes 

No closed-form inversion of (22,23) from D, to the primitive variables Po,u, ft* is known, so 
the primitive variables must be found numerically. In this section, we describe and discuss several 
methods for solving these equations. 

A popular and robust means of solving systems of algebraic equations numerically is the 
Newton-Raphson (NR) scheme (see Section 9.6 of Press et al. (1992)). Since it is also simple to 
code, we will employ this method by default. 


3.1. IDw Scheme 


It seems likely that the solution would be obtained more efficiently if one could manipulate 
the five fundamental equations to reduce the dimensionality of the system, as was done by Del 
Zanna, Bucciantini, & Londrillo (2003) for special relativistic MHD. Our procedure for reducing 
the 5D system is to evaluate certain scalars from the conserved variables, then expand expressions 
for these scalars to simplify the solution. We will use D, and which is the energy 

density measured in the normal observer frame. We will also use where Q'^ = j ’^which 
is the energy-momentum flux perpendicular to the normal observer. Ultimately, we obtain two 
equations dependent only on the known conserved variables and two independent variables 7 and 
w. Instead of these two unknowns, however, we use W = and = UjU*, where u* = u^/'j is 
the flow velocity relative to the normal observer. These new variables simplify the equations and 
lead to numerical schemes that are more robust. Whenever 7 is stated from now on, it is implied 
that it is calculated from via the identity 7 ^ = 1 / (l — u^). 


A consistent procedure can be developed as follows. First, expand the definition of Qf^ to find 
the key result 

B^^Q^ = Wh. (26) 

Since B^Qfj, can be evaluated in terms of the (known) conserved variables, this equation can be 
used to eliminate Ufj_B^ in favor of the unknowns and W. 


Next, expand using (26) to find 


7 = + 


Finally, we solve for to find the simple result 


u2 = 


Q 2^2 ^ (^2 ^ 2IU) 

{B^ + Wf IU2 


(27) 


(28) 


This is the relativistic counterpart of the nonrelativistic expression where P is the 

magnitude of the momentum density vector pv. It gives us an explicit expression for v'^iW) that 
results in positive definite evaluations if implemented in the code as shown here. 
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Next, expand Qfj,n^ using (26) and the variables 7 and W: 

Q ij,n^ = — — {1 + v‘^) + W+p{u,po). (29) 

This remaining equation then yields a solution for W. Specifically, the 1D\y scheme solves one 
nonlinear algebraic equation (29), which is now only a function of W since (28) is used to eliminate 


Once and W are found one can recover w, po (from D), and u. The next step is to find u 
using Q*. After some manipulation one finds 


Q^^ = -{W + ^2) u^^ - 

7 7 

Since u^B^ = {'y/W){Qfj,B^), this can be used to solve for ft*: 


u = 


7 


W + B^ 




w 


(30) 


(31) 


3.2. 2D Scheme 

The equation to be solved in the ID^v method includes a quotient of polynomials in W since 
it implicitly uses (28) for v^. This suggests that numerical pathologies might arise near roots. By 
solving the two, simpler equations (27,29) simultaneously for W and v^, one may eliminate such 
problems. We call this method the 2B scheme since it involves solving a two-dimensional algebraic 
system. 

We find that using instead of or 7 is particularly advantageous for this method. This is 
because equations (27,29) are linear only in (modulo the u^-dependence of the state equation) 
and not in or 7 . The linear dependence on increases the rate of convergence for this quantity 
and is guaranteed to be well-behaved in the vicinity of a root. 

Koide et al. (1996) and Koide, Shibata, & Kudoh (1999) also used a two-dimensional method. 
But instead of and W, they use (7 — 1) and {u^B^) as independent variables. We have not 
tried their method since the functions they minimize are more complicated than ours and assume a 
T-law state equation. Further, it is likely that these two methods would perform similarly since one 
can eliminate W for u^B^ in our method via equation (26). As mentioned earlier, we find better 
performance using instead of 7 . 


3.3. 5D Scheme 

The simplest procedure, and the one we used initially in the HARM code, is to invert the five 
equations (22) to (23) using a multidimensional NR scheme. This requires evaluating the matrix of 



derivatives 3U/(9P. Further, a Newton iteration of this system involves more operations than the 
IDw and 2D schemes since it requires calculating elements of d\J/dP and a matrix inversion. Also, 
we find that it requires an initial guess that is close to the solution. This is almost always available 
from the last timestep (and if the guess is not good it usually means something is wrong). Because 
it involves root-finding in a five-dimensional space we call this the 5D scheme. The conserved 
variables used in this method are the same as those used in Gammie, McKinney, & Toth (2003). 
All other methods described here use D and Q^. 

Notice that the ID^y, 2D, and 5D schemes do not require a particular equation of state. 
For example, derivatives of p with respect to the independent variables could be obtained from 
equation of state tables using finite difference approximations. Two of the next three methods, 
however, assume a F-law equation of state. 


3.4. ID .^,2 and 1D*2 Schemes 


Another scheme can be derived if we assume that the equation of state is 

/ ^ r-1/ D\ 

p = (r-l)u = j. (32) 

If we make this substitution and use the expression for v^{W), equation (29) reduces to an eighth- 
order polynomial in IT, for which there is no general closed-form solution according to the theorem 
of Abel-Ruffini. It is simplest to solve this single nonlinear equation numerically. Because of the 
complexity introduced by the equation of state there is likely no general way of isolating the physical 
root, and one must simply look for a solution that is close to the solution of the last timestep. With 
other state equations, it may not be possible to express (29) in polynomial form at all. 


It is worth noting how this situation resolves itself in relativistic /lydrodynamics. Equation 
(28) becomes 


u^ = 


and equation (29) becomes 


IT2 - Q2 


=p{W,Dli)-W. 


(33) 

(34) 


Obviously, there is no general solution since one has the freedom of choosing an equation of state. 
Our particular equation of state, however, yields a quartic whose solution was discussed by Eul- 
derink & Mellema (1995). 


Instead of solving the eighth-order polynomial directly, the scheme, which is a modified 

version of the special relativistic method described in Del Zanna, Bucciantini, & Londrillo (2003), 
solves a cubic equation for W{u^) and a nonlinear equation for u^. The cubic described in Del 
Zanna, Bucciantini, &: Londrillo (2003), however, can sometimes have two positive, real solutions 
for IT. The larger root appears to be the physical one, though a general proof of this has eluded us. 
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In order to eliminate this uncertainty, we use the following cubic equation which we have proved 
leads to only one, positive solution^ 

C{W) = W^+ 3d2W^-Mo = 0 (35) 


where 


8[i + u 2 (r- 1 )] 


(36) 


d 2 = 3 ^2 (P _ ]_)j + -^ (l + + -C* (1 - 1/r) \/l - . (37) 

Equation (35) results from multiplying equation (29) by FW'^/ [l + (E — 1)] and substitut¬ 
ing equation (32) for p(W, D/j). In general, there are three solutions which are given in closed-form 
by Cardano’s formula (Weisstein 2004; Rade and Westergren 1990). If our cubic has a positive and 
real solution then there is only one positive and real solution, and it can be shown that this solution 
is always equal to the following 



where 

V = Ado {do - dl) . (39) 

It is useful to know that do > 0 always and d 2 can take negative and positive values. When V < 0, 
the solution can be expressed in a simpler form: 

W{V < 0) = ds [cos (d/3) - 1] , 9 = cos-^ [2do/d^ - l] . (40) 

Table 1 lists the possible physical solutions of (35) depending on the particular values of do and 
d 2 - In the 1D„2 scheme, we use the physical solution for W{v‘^) and numerically solve an equation 
proportional to equation (27) for That is, 1D^2 solves for via a NR method in which the 
physical cubic solution is calculated for each iteration. 

One drawback of this method is that it requires that the state equation be linear in IT. Solving 
equation (29) numerically instead makes the method compatible with general equations of state. 
We call this technique the 1D*2 scheme. It consists of taking NR iterations to find IT(u^) nested 
within Newton iterations to find v^. In other words, for each Newton update of we solve equation 
(29) for IT using a separate NR method assuming the most current value for v^. This supplies 
the next iteration with a consistent value for IT. Surprisingly, this nested NR method (1D*2) 
is faster, more robust and more accurate than the 1D^2 method. The difference in accuracy and 
efficiency is likely due to the appearance of transcendental functions and condition statements in 
the closed-form solution of (35). 


^The signs of the cubic solutions can be derived from the locations of the cubic’s local extrema. The existence of 
only one positive solution stems from the following properties of the cubic: 1) dC/dW\w=o = 0; 2) C{W = 0) < 0 
always; and 3) limw^ioo C(W) = Too. 
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3.5. Polynomial Scheme 

The substitution of equation (28) into equation (29) leads to an eighth-order polynomial in 
W if we assume a T-law equation of state (32). This suggests the possibility of using a general 
polynomial root-finding method (such as Numerical Recipes’ zroots) that finds all 8 roots. We 
will call this the polynomial scheme. 

The physical root can be identified by requiring that it also solve the five equations U = U(P). 
Unfortunately this test can sometimes yield ambiguous results due to amplification of roundoff error, 
making it difficult to identify the correct solution. This method also turns out to be computationally 
expensive. 


4. Numerical Implementations 

In the previous section, we described the mathematical framework that embodies the primitive 
variable inversion methods we have tested. We will now discuss the details of the numerical methods 
we have used to test the performance of these formulations. All routines (in the C language) are 
available for download from the web.^ 

We use the Newton-Raphson scheme to solve the nonlinear algebraic equations of the IDw, 
2D, 5D, 1D*2, and 1D„2 methods. An excellent description of this method—as well as other root¬ 
finding methods—can be found in Press et al. (1992); we will only state unique aspects of our 
implementation here and defer further explanation to our source code. Let R denote the system 
of nonlinear equations, or “residuals,” x the independent variables for which we are solving, and J 
the Jacobian whose (i, j)-component is dKi/dxj. The NR procedure assumes that R is smooth and 
nearly linear, so that we can make consecutive linear corrections to our guess that lead us toward 
the root. This NR step is defined as the following; 

Ax = -J"^ • R . (41) 

We measure convergence using a “Newton-Raphson error function,” Snr = |A1U/1U|, where 
AW is the change in W between the two most recent NR iterations. We find that other convergence 
criteria yield little if any improvement, which is not surprising since W is the crucial variable. In 
particular, a convergence criterion based on residual errors in the solution of equations (22,23) is 
not used for any of the schemes since it is difficult to normalize the error properly a priori over the 
entire parameter space. 


^The current version of the source code for each method described in this paper can be 
found in the electronic edition of the Journal, while a maintained version can be obtained at 
http://rainman.astro.uiuc.edu/codelib/codes/pvs_grmhd/. The methods can be used with any user-supplied 
spacetime metric. 
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The condition for stopping the iterative procedure involves three parameters: TDL, Ntsih and 
-^extra- If -®NR falls below TOL after performing at least one iteration, then only A^extra more 
iterations are performed. A solution is said to be found if the tolerance criterion is satisfied before 
IVnr iterations. 

The impetus for the additional Aextra iterations is one of efficiency and accuracy. When the 
error reaches our tolerance level, the extra iterations try to reduce it even further. Since this can 
sometimes be fruitless—e.g. the solution error becomes insensitive to subsequent NR steps—we 
set Aextra to a Small number. The parameters used by default in this paper are TDL = 10“^^ and 
Aextra = 2, which should yield solutions accurate to within roundoff error if the guess is sufficiently 
near the true solution, since the NR method converges quadratically. 

What value should be chosen for TDL? Ideally the inversion error TDL should be smaller than 
the truncation error that is inevitably introduced in the numerical evolution of U. But truncation 
error is difficult to estimate on-the-fly ^ and cannot usually be estimated a priori. 

Different problems may require different values of TDL. One class of problems that is particu¬ 
larly sensitive to the value of TDL is the evolution of small amplitude waves. This is a problem of 
some interest, as it is frequently used in convergence testing numerical methods. In this case the 
fluctuating part of the fluid variables is small, and the truncation error is a small fraction of that. 
The truncation error can therefore be a very small fraction of the conserved variables, and TDL 
must be comparably small to avoid spoiling convergence. On the other hand, the measurements of 
the accretion rates of rest mass, energy, and angular momentum in the relativistic disk simulation 
discussed below appear to be peculiarly insensitive to TDL. The only safe course is to directly check 
the sensitivity of numerical measurements to TDL. 

During a sequence of NR iterations, an independent variable may leave its allowed domain 
(e.g. superluminal velocities). In order to prevent numerical divergences and unwanted imaginary 
parts, we reset the independent variable to a value in its physical domain. For example, in the 
IDw, 1D^2, 1D*2, 2D and polynomial methods, we demand that 0 < < 1 and IT > 0. No 

constraints are used in the 5D method, since experimentation with this has led us to believe that 
constraining po,u > 0 only increases the likelihood of not converging to a solution. 

A special note on the polynomial scheme, where we use Laguerre’s method to find IT: there 
is an additional numerical criterion that must be used for evaluating which of the 8 roots is the 
physical root. One first eliminates roots that have imaginary components above some threshold 
value, then evaluates the residuals in the solution of the basic equations for the remaining roots. 
The root with the smallest residuals is identified as the physical root. As for the NR methods, 


®Local truncation error estimation can be accomplished in place during a simulation without knowing the exact 
solution beforehand by solving the equations of motion additionally on a “shadow” grid whose grid spacing is half 
that of the base grid. The truncation error estimate is then calculated by comparing the two solutions at coincident 
timesteps (Berger & Oliger 1984; Choptuik 1995). 
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parameters TOL, A^nr and A^extra are used to control the accuracy of the solutions. 

A few final comments on numerical implementation. In terms of complexity, the 1 D* 2 , 

2D, and polynomial methods are the easiest to implement; some care is required in finding the 
physical solution to the cubic equation in the 1 D ^2 scheme. The lower-dimensional NR schemes 
are simple enough to derive the closed-form expressions for R, J and Ax by hand, so we have 
coded those in directly. The 5D scheme, by contrast, is complicated by the need to evaluate the 25 
elements of J. This can be done analytically (this procedure is time-consuming and prone to error) 
or numerically using finite differences (this is simple but introduces additional numerical noise). 


5. Tests 

The space of possible numerical approaches to the inversion problem is large, so we can make 
no claim that any of our methods are optimal or near-optimal. In this section, through a series of 
tests, we show that some of the methods are “good enough” in the sense that (I) they can be shown 
not to contribute significantly to the numerical error in actual astrophysical applications, and ( 2 ) 
some of them are efficient enough that they do not contribute significantly to the computational 
cost of an evolution. We have expended significant effort in optimizing each scheme’s speed and 
accuracy, so that these tests make a fair comparison between the various methods. 

We consider three tests. The first is a parameter space survey in which the primitive variables 
are varied over many orders of magnitude and the accuracy of the solution is evaluated. The second 
test places the inversion routine in a special relativistic MHD code and evolves a magnetized, cylin¬ 
drical explosion problem due to Komissarov. The third test places the inversion routine in a general 
relativistic MHD code and evolves a magnetized disk around a rotating black hole. Throughout we 
assume a the T-law equation of state (32) with T = 4/3. 


5.1. Parameter Space Survey 


The inversion routine takes as arguments the 8 conserved variables (three are magnetic field 
components and are trivially converted to primitive variables), guesses for the 5 nontrivial prim¬ 
itive variables, and the 10 components of the metric. This parameter space is too large to cover 
exhaustively, so we only vary po, u, 7 (the Lorentz factor), and = cos“^ (^iB^/^UiU^BjB^^ 
over the region spanned by a typical accretion disk evolution; 


logio Po G [-7,1] > logio^x G [-10,0] , logio 7 G [0.002,2.9] > 
logio ^ [-8,1] , cos 4> G [-1,1]. 


(42) 


We sample uniformly in the logarithm for each variable. Specifically, 40 x 40 x 20 x 20 x 9 points 
are evenly sampled along dimensions log^g /^o < 8 ) logio ^ ® logio 7 logio ® respectively. In 
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order to choose reasonable relative magnitudes between the components of h* and at a given 
location with respect to the Kerr-Schild metric, we select 9 points from an accretion disk simulation 
(see Section 5.3) at which cos(^>) ~ —1,..., 1. The specific values of cos(<h), ti*, and coordinates 
used for the parameter space survey are given in Table 2. The overall magnitudes of the u* and B^ 
are set by the values of 7 and B^ at the given parameter space point. 

At each point of the parameter space we perform the forward transformation; this gives a 
(nearly) exact solution to the inversion. We then feed the inverter a guess obtained by multiplying 
each primitive variable by 1 + d, where d is a random value between —1 and 1. The same sequence 
of random values is used for each method. This turns out to be quite a stringent test, particularly 
when the random value is near —1. Even though the maximum threshold for the offset may not be 
large enough to model the behavior at strong shocks, it is approximately the maximum seen for po 
and u in the bulk flow of our accretion disk simulations (see Section 5.3). 

Our parameter space surveys were made using TOL = Ntsir = 30 and Ng^tra = 2; these 

values were determined after the fact to be nearly optimal. 

Figure 1 shows the normalized error in the internal energy, Au/u, for all the methods. These 
accuracy measurements indicate how strongly roundoff errors are amplified by the scheme. The 
error is averaged over all the parameters except u. Only the points for which all the methods 
converge are included in the average^. Two methods stand out as less accurate: the polynomial 
scheme, and the 1 D ^2 scheme. An examination of Au/u over the entire parameter space indicates 
that the the IDw, 2D and 1 D *2 methods yield almost indistinguishably accurate solutions. If one 
ignores failures points, the 5D method is often as accurate as these methods. The 1 D „2 is the fifth 
most accurate (ignoring failure points), and the polynomial scheme is the least accurate. 

The relative errors in po and id increase with increasing B'^/po, B‘^/u, and 7 for the most 
accurate schemes. The relative error of u —for the same methods—grows with po/u, 7 and B'^/u, 
and is fairly independent of B'^/po- We find that there is typically a maximum value of po/u above 
which the methods nearly always fail; this maximum value starts large (~ 10 ^^) but reaches unity 
by the time 7 ~ 10® We find that poor accuracy (> 1% relative error) in po and u* usually 
occurs for points with B'^/pa, B‘^/u > 10^®, but u is accurate for B'^/u < 10^^. The dependency on 
7 seems to not be strongly tied to any of the other variables, and we find that—on average—> 1 % 
relative errors are seen when 7 > 10 ^ — 10 ®. Also, the accuracy of all the methods improves as 


^Since the 5D method fails at almost half the original points and a majority of these failures occur when 7 is 
large, this may introduce a bias. We have performed an identical survey in which no points were neglected. Even 
though solution rates decrease slightly, the rankings in Table 3 stay the same. The high-ENR tails (Figure 2) extend 
to larger values of the Enr distributions of the 2D, IDw, 1 D() 2 , 1 D „2 and 5D methods now extend to ~ 

10“®, 10“®, > 1 and > 1, respectively. More iterations are required and all methods but the 2D method have A^nr 
distributions extending to ATtr = 30. The 2D method remains the best method. 

®These numbers were obtained in an extended parameter space survey in which we sampled 40^ points for each 
of the 9 values of $ over the space logjgfpo), logiQ(rt), log^gfi?^) G [—10,10] , and logj^pfy) G [0,6]. 
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cos <I> ^ 0, on average. Note that the precise values of these thresholds are not universal to all 
situations, and are quantitatively dependent on all the degrees of freedom, including the offset d. 

The inversion routine is said to “fail” if it has not converged after iV^R = 30 iterations. Failure 
rates for the methods are given in Table 2. Here the 2D method stands out as the best, failing only 
5 times in the survey of 5.76 x 10® points. The VDw and 1D*2 methods also have low failures rates. 
The 1D^2 and polynomial schemes fail much more frequently, at a rate about an order of magnitude 
lower than that of the least robust 5D method. Notice that while the polynomial method converges 
more often than the 5D method, it often converges to an inaccurate result. 

We have found that all the methods asymptote to a minimum failure rate as the range in 
d is diminished, i.e. as (5i ^ 0 when d G [—5i,5i]. The 5D method’s asymptotic regime ends 
at approximately (5i ~ 10“^, while the other methods are insensitive to changes for (5i < 1, where 
= 1 is the maximum value for which the po and u guesses remain non-negative. If we instead vary 
the overall magnitude of the guess, i.e. (1 -|- d) ^ 82 (1 -|- d), set di = 1, and set 82 G [10“®, 10®], 
we find that the 5D method has a minimum failure rate at about ^2 ~ 1 as expected. The other 
methods fail at different, yet constant, minimum rates up until 82 ^ 100; this may be because 
the residuals used in these routines steepen as VF ^ 00 . This does not imply that every time the 
density jumps by ~ 10^ the method will fail, since the relative offsets are not likely to be large for 
all primitive variables at the same time in practice. The high failure rate and greater sensitivity to 
the guess’ offset seen with the 5D method is easily attributed to the fact that higher-dimensional 
minimization problems are more difficult than lower-dimensional problems (e.g. Press et al. (1992)); 
more dimensions admit residuals with more complicated landscapes and make a method more likely 
to be mired by local minima. 

As another means of demonstrating the accuracy of each method’s solutions, we show his¬ 
tograms in Figure 2 of the parameter space points binned against the method’s exit value of F^nr- 
Any points for which Anr < 10“^® are considered acceptable solutions; since only those points for 
which solutions are found by all methods are included, there are no points beyond F^nr = 10“^®. 
Those methods that explicitly solve for W (i.e. 2D and IDw) give rise to a power-law distribution 
extending up to machine precision since the change in W between iterations AW is explicitly cal¬ 
culated during the NR step. The other methods result in non-zero AW/W only when it is above 
machine precision since they calculate it by subtracting the previous iteration’s value from the 
present one. From a numerical perspective, anything below machine precision should be considered 
equivalent to zero anyway, so we shall only concern ourselves with the distributions above machine 
precision. In this regime, the 2D method again yields the best results, followed by 1D*2, IDr/, 
1D^2 and 5D in order of increasing average F^nr. 

What is the maximum accuracy attainable by each method? We have surveyed parameter 
space using values of TDL = 10“^® and Nnr = 200, which should give close to maximum accuracy. 
The distributions of F^nr for the IDyi/, 2D and 1D*2 methods remained nearly unchanged from 
those shown in Figure 2 which used the default (less accurate) set of parameters; the additional 



- 15 - 


iterations failed to significantly reduce solution error, on average, for these methods. The 5D 
method, however, does benefit from a larger iVNR. Its accuracy and failure rate improve the longer 
it is allowed to iterate, although it still underperforms the 2D, IDr/ and 1D*2 solvers. 

In order to measure how quickly the methods converge to a solution, we have plotted histograms 
in Figure 3 of the number of NR iterations, Nnr, each method performed before exiting per 
parameter space point. From the figure, we see that the ID^v, 2D and 1D*2 methods all typically 
use between 6 — 10 iterations, and that they give rise to nearly Gaussian distributions in this 
regime. The 1D^2 method has an extended power-law tail, while the 5D method has a nearly 
constant distribution over Nnr. Even though the 1D^2 scheme has a peak at small Nnr that is 
narrower and taller than the other schemes—suggesting that it may be more efficient— it fails to 
find a solution quite frequently and has a substantial tail. The IDr/ method appears to have the 
best convergence rate over the sampled parameter space since it has a steep exponential tail and it 
usually converges to a solution in less than 10 iterations. 

Which scheme is fastest? Table 3 shows the number of solutions per second achieved by each 
scheme in the parameter survey, on our 3.06 GHz Intel box. The two fastest methods—ID^v and 
2D—are about 30% faster than 1D^2 method, which also has the smallest average Nnr. Since the 
1D*2 method performs many more computations (because of its nested NR scheme), it takes fourth 
place. The 5D method is almost an order of magnitude slower, due to its larger linear system and 
larger average Nnr. The polynomial method is a factor of 2 slower than the 5D. 

We have tried a variety of schemes for improving the accuracy and convergence rate of each 
scheme. For the NR schemes, we tried implementing a line searching method, which attempts to 
find an optimal step size along the direction of steepest descent (Press et al. 1992). This method 
saves a few of the nonconvergent outliers in the parameter space survey for the IDr/ method but 
does not help anywhere else for any other method. We tried this improved IDr/ scheme in an 
accretion disk simulation and found that it did not improve the success rate at all in the disk 
simulation. 

To sum up, the 2D method has the lowest failure rate and is only slightly more computationally 
expensive than the fastest method. If the initial guesses are close enough to the solution, the 5D and 
1D„2 methods fare nearly as well as the more successful IDw , 2D, and 1D*2 methods in accuracy 
but because of their larger rates of nonconvergence neither is recommended. The polynomial 
method is slow and inaccurate, and we do not recommend it. We will ignore the worst of these 
methods in the next two sections and discuss the performance of only four methods; 2D, IDr/, 
1D*2, and 5D. 


5.2. Cylindrical Explosion 

The parameter space survey may emphasize different aspects of the inversion routine than a 
realistic RMHD problem. Here we consider inversion routine performance in a special relativistic, 
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slab-symmetric, magnetized explosion. We use the same initial conditions as Komissarov (1999). 
The computational domain is defined in Cartesian coordinates (rc, y) G [—6, 6] x [—6, 6] discretized 
by 200 X 200 points. A pressure-enhanced cylinder is centered on the origin with radius r = 0.8, 
and is matched exponentially to a constant background that starts at r = 1. The inner state has 
Po{r < 0.8) = 10“^ and P{r < 0.8) = 1, while the outer state is initialized with po{r > 1) = 10“^ 
and P{r > 1) = 3 x 10“®. All spatial components of the velocity are zero at t = 0, and the magnetic 
field is uniform, pointing in the x-direction: = (0,0.1,0,0). 

We evolved the initial state using a version of the HARM (Gammie, McKinney, & Toth 2003) 
code. We use an HLL (Harten et al. 1983) flux, and linear (second order) reconstruction with a 
minmod slope limiter. We use a Courant factor of 0.1 in the evolution. Figure 4 is a reconstruction 
of Figure 10 of Komissarov (1999) using our own evolution; the results are qualitatively identical. 

All the inversion schemes perform similarly on this problem. Table 3 outlines the speed and 
failure rates obtained with each method. When a solution is found, it is found by all methods 
in less than 10 iterations. The only significant differences are in the run time: the 2D is clearly 
fastest, while 5D is a factor of 2 slower than the next fastest method. The failure rates are the 
same, because the methods fail at the same points where the U calculated from the evolution is 
unphysical. The small Courant factor may play a role in the low failure rate here, since smaller 
timesteps imply smaller changes in the U over a timestep, so the initial guess used in the inversion 
routine (P from the previous timestep) is closer to the solution. 


5.3. Accretion Disk Evolution 

A more challenging context for testing the inversion routines is the evolution of a magnetized, 
centrifugally supported torus around a Kerr black hole (see, e.g., McKinney & Gammie (2004) and 
Gammie, Shapiro, & McKinney (2004) for similar applications). In this particular incarnation of 
the problem we evolve a Fishbone & Moncrief (1976) torus containing a weak poloidal seed field 
around a Kerr black hole with a/M = 0.9375. The initial state is unstable to the magnetorotational 
instability (Balbus & Hawley 1991), so turbulence develops in the disk and material accretes onto 
the black hole. 

The inner edge of the initial torus lies as QGM/c^ and the pressure maximum lies at VIGMj(?. 
The Fishbone-Moncrief equilibrium condition of u^v!' = 4.28 is used. On top of the Fishbone- 
Moncrief torus we add a weak magnetic field with vector potential = Max (po/pmax — 0.2,0) 
where pmax is the maximum of the disk’s rest-mass density. The magnetic field amplitude is 
normalized so that ratio of gas to magnetic pressure within the disk has a minimum of 100. With 
the addition of the field the disk is no longer strictly in equilibrium, but it is only weakly perturbed. 

The initial conditions are evolved using the HARM code. The flux is a local Lax-Friedrichs 
flux, and linear (second order) interpolation is used with a Monotonized Central limiter. Since 
HARM is incapable of evolving a perfect vacuum, we surround the disk in an artificial atmosphere. 
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or “floor” state, with /3o,atm = 10 ^{r/M) and Uatm = 10 ^{r/M) Also, po and u are set 
to their floor values if and when they fall below the floor at any point in the evolution. 

The equations of motion are solved using finite difference approximations on a discrete do¬ 
main of so-called Modified Kerr-Schild (MKS) coordinates. Instead of the standard Kerr-Schild 
coordinates {t,r,0,(p), a uniform discretization of MKS coordinates is used 

where 

r = , 6 = + ^{1 — h) sin(27rx*'^^) . (43) 

This coordinate transformation is intended to concentrate cells near the equator and at small radius. 
The cells are centered at coordinates 

xf ^ ^ 

where i € [0, A^i — 1], j G [0, A "2 — 1], and 

= log (nn) , ^ log f—] , ^ . (45) 

JVi \ Tin / JV2 

A^i and N 2 are the number of cells along the x^^l and x^^l axes, respectively; for the 128^ (256^) 
run, Ni = N 2 = 128 (Ni = N 2 = 256). rjn is chosen so that approximately ten cells lie behind the 
event horizon: rjn ~ 1.0035M for the 128^ run, and rjn ~ 1.1702M for the 256^ run. The remaining 
parameters are always set to rout = 40M and h = 0.3. 

We have run this torus problem using the 2D, ID^v, 10^2, and 5D inversion schemes. The 
parameters are the same as before: {TDL, Anr, Aextra} = {10 ^^,30,2}. The resolution is 256^. 
The initial conditions in each run also contain identical low level noise to initiate the growth of the 
magnetorotational instability. The accretion rates of rest-mass, energy, and angular momentum 
for these runs are shown in Figure 5. One would expect that after the disk becomes turbulent 
(at ~ 500 GM/c^ in this model) small differences in the solution due to the inversion routines 
would be amplified and that the accretion rates would not track each other very well. Evidently 
the accretion rates are nearly identical until 900 GM /, after which their rates follow only vaguely 
similar trends. This suggests that the differences in evolution due to the inversion routines were 
small indeed. 

The average normalized accretion rates for the internal energy (E/Mo), angular momentum 
(L/Mo) and rest-mass (Mo) given in Table 5 differ by only a few percent. Far larger differences 
were obtained by changing the seed used to generate the noise in the initial conditions. In addition, 
the qualitative structure of the disks in steady-state remains the same at the end of the evolutions. 
Shown in Figure 6 are snapshots of log(po) at t = 2000 GM/(? from the runs using the 2D, IDw, 
1D*2 and 5D methods (shown left to right, respectively). The most pronounced differences are in 
the corona and funnel regions overlying the bulk of the disk. Evidently the solution is stable to 
changes in the inversion algorithm. 
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We have also compared low resolution (128^) runs using the four methods. This is interesting 
because P exhibits larger fractional changes per timestep at lower resolution, so the inversion 
routines are put under greater stress. In Figure 7, we show the distribution of i?NR exit values over 
the course of these simulations. The most striking feature of the plot is in the number of points 
at which the 5D method fails to converge compared to other methods. The 2D method rarely fails 
to converge, and the IDw and 1 D *2 methods fail to converge at nearly equal, intermediate rates. 
In Table 6 we state the failure rate—i.e. the frequency at which either an unphysical P value is 
found or no solution is found—for each 256^ run. The table demonstrates that the failure rates 
are approximately the same, suggesting that when the 2D method “fails” it is converging to an 
unphysical P solution instead of failing to converge altogether like the 5D method. 

The inversion schemes are also faster in the disk evolution than in the parameter space survey. 
This is seen in the distribution of A^nr over the lower-resolution run shown in Figure 8. Many 
more inversions in the disk run converge sooner than in the parameter space survey. The greater 
efficiency and accuracy seen in the disk evolution is likely due to the fact that in the disk evolution 
the inversion routine is usually supplied with better initial guesses, from the previous timestep, 
than in the parameter space survey. 

What is the optimal value of the Newton-Raphson parameters, {TOL, Nnr, iVextra}, for the 
most successful 2D and methods? Almost no new acceptable solutions are found if we increase 
.^NR > 30, so this is the natural value for this parameter. We also reran the disk evolution using 
TOL = {I0“^, 10“^, 10“'^, I0“®, I0“^°, 10“^^} at a resolution of 256^. Each run was otherwise 
identical and each used {Anr, Aextra} = {30,2}. Surprisingly, we discovered that the evolutions 
deviated little from each other until after the inner parts of the disks become turbulent, at f ~ 
500M. The relative difference of any given primitive function between any run and the TOL = 10“^^ 
run grows as ~ until it plateaus to a constant value at f > lOOOM. Further, the accretion rates— 
i.e. those functions displayed in Figure 5—were qualitatively indistinguishable until t ~ 800M. The 
total number of NR iterations executed during the course of a simulation also varied little between 
these runs. 

Similar runs were made for TOL = {I0“^, I0“®, 10“^*^, 10“^^}, but now we reduced the reso¬ 
lution to 128^. Since the timestep is larger in the lower resolution runs, guesses for the primitive 
variables are—on average—further away from their solutions than in the 256^ resolution runs. As 
a result, the accretion rates and distribution of Enr differed dramatically between the TOL = 10“^ 
run and the others. Runs with 10“^^ < TOL < I0“® have similar accretion rates, and distribu¬ 
tions of Enr and Anr. This implies that the guesses in the higher resolution run are so close to 
their solutions that performing 3 iterations—the minimum number of iterations allowed in the NR 
schemes when Aextra = 2—often leads directly to an accurate solution independent of TOL. 

These tests suggest that the tolerance, if below at least 10“^, is inconsequential to our accretion 
disk simulations when the grid resolution is sufficiently high and when Aextra = 2. To evaluate the 
sensitivity of the outcome to Aextraj we ran two 128^ disk models with TOL = 10“®; one with 
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no extra iterations and the other with A^extra = 2. The quantities L/Mo, E/Mq, and Mo deviate 
between the two runs by less 4% at any given time, far less than what is seen between two otherwise 
identical runs at resolutions of 128^ and 256^. We conclude that the inversion error is in some mean 
sense much less than the discretization error when TDL is below 10“®. 

The -Enr distribution of the A^extra = 0 run is, however, significantly different than the one 
shown in Figure 7. Without the extra iterations, the -Enr distribution is approximately uniform 
over 10“^® < Enr < 10“®, and becomes more similar to the A^nr = 2 distribution for Enr < 10“®. 
The extra iterations therefore seem to be doing what we expected for the most part. In addition, 
using two extra iterations per inversion only increases the simulation’s run-time by 4%. Even though 
this doubles the run-time contribution from the primitive variable solver, it is still an insignificant 
portion of the evolution’s total computational cost, which seems worthwhile to nearly eliminate the 
possibility that inversion error makes a significant contribution to the computational error budget. 

This investigation highlights the importance of handling aberrant cells in real applications. 
When an inversion fails to converge to a solution, when a solution leads to unphysical primitive 
variables, or when it converges to a set of primitive variables which we have some a priori reason for 
classifying as numerical artifacts, one must either halt the run or else “correct” these values®. For 
example, the most common problem involves evolving to a state of negative internal energy density. 
We have found that using an artificial floor can lead to the generation of spontaneous “explosions.” 
These explosions can corrupt and eventually terminate the evolution. We have found that 2“^'^- 
order interpolation of P from the problematic cell’s nearest-neighbors is successful at eliminating 
these cell-scale artifacts. This interpolation procedure is used whenever at least one of the above 
conditions is met. All disk evolutions mentioned in this paper used this method. 


6. Conclusion 

We have outlined a compact derivation of the equations for calculating P(U). This formulation 
suggests a variety of possibilities for performing the calculation numerically. We have implemented 
a small subset of the possible methods, and compared these implementations through (1) a survey 
over a large subset of possible U values; (2) embedding the inverter in a SRMHD code and evolving 
a cylindrical explosion problem due to Komissarov; (3) embedding the inverter in a GRMHD code 
and evolving a turbulent, magnetized disk around a rotating black hole. 

Several key points emerge from the comparison. First, the implementation can be made ac¬ 
curate enough that variable inversion does not make a significant contribution to the code error 
budget. Second, some implementations can be made fast enough that they occupy only a few 
percent of the total cycles used by our GRMHD code. Since we are using a particularly sim¬ 
ple GRMHD algorithm, this is likely an upper limit to the fractional cost of the inversion in all 


®In accretion disk simulations we typically classify a point as aberrant if po,M < 0, 7 > 50, or 7 < 1. 
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GRMHD schemes. Third, the inversion routine originally used in HARM (Gammie, McKinney, & 
Toth 2003), called “5D” here, is particularly slow and inaccurate. 

We recommend the “2D” scheme, described in §3.2. A version of this scheme is available for 
download on the web at http://rainman.astro.uiuc.edu/codelib/codes/pvs.tgz 

This work was supported by NSF grants PHY 02-05155, AST 00-93091 and PHY99-07949 
(Kavli Institute for Theoretical Physics, where a portion of this work was completed). We thank 
Yuk Tung Liu, Po Kin Leung and Xiaoyue Guan for comments. 
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Table 1. Physical Solutions for W of the 1D^2 Scheme’s Cubic 


Case 

Condition 

V 


1 

dl> do >0 

P < 0 

(40) 

2 

dl < do > 0 , (i2 / 0 

P > 0 

(38) 

3 

d2 = 0 , do ^ 0 

V>Q 

W = (4do)^/^ 

4 

d2 = do > 0 

T> = 0 

II 

5 

do = 0 , ^2 / 0 

T> = 0 

IT = -3d2 iff d2 < 0 

6 

to 

II 

o 

II 

O 

T> = 0 

(none)'’ 


'^Only physically-acceptable solutions, i.e. those that are real 
and positive, are given here. 

'’'The only solution, IT = 0, is unphysical. 
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Table 2. Kerr-Schild Coordinates (a = 0.9375) used in the Parameter Space Survey 


cos(^>)^ 

r 

0 

-0.751 

8.195 

1.552 

-0.250 

1.375 

1.444 

-0.500 

2.676 

1.016 

1.000 

23.166 

2.672 

-0.997 

26.467 

0.658 

0.500 

1.571 

1.589 

0.749 

3.588 

1.455 

0.250 

2.406 

2.483 

-0.0005 

35.480 

0.146 


Note. — Please refer 
to the electronic version 
of the Journal for a com¬ 
plete, machine-readable 
table of the parameters 
used for the survey. The 
larger table provides 

for each value of cos <1>. 
Please refer to Section 5.3 
for a definition of the 
coordinates. 


®'cos(<h) 


UjB'- 
^UiU'^ Bj 



Table 3. Parameter Space Efficiency Comparison 


Method 

NR steps per sol. 

Sol. per sec. 

Failure Rate 

2D 

8.45 

1.66 X 10® 

8.7 X lO-"^ 

IDw 

7.45 

1.68 X 10® 

8.8 X 10“^ 


7.08 

1.06 X 10® 

3.6 X 10“^ 

1D^2 

7.05 

1.24 X 10® 

2.0 X 10“2 

5D 

19.3 

1.89 X 10^ 

4.2 X 10“® 

Poly 

— 

9.21 X 10® 

4.1 X 10“2 


Note. — Average number of NR steps taken per solution, 
average solution rate, and average failure rate (per solution) 
are shown for each method. The first entry for the polyno¬ 
mial method is vacant since it does not use the NR scheme. 
The Intel C Compiler Version 8 and an Intel Xeon 3.06GHz 
workstation were used for these runs. 


Table 4. Cylindrical Explosion Efficiency Comparison 


Method 

NR steps per sol. 

Zone-cycles/sec. 

Failure Rate 

2D 

3.81 

74142 

3.75 X 10-'^ 

IDw 

3.76 

68966 

3.75 X 10"^ 


4.68 

58042 

3.75 X 10"^ 

5D 

4.45 

27100 

3.75 X 10-^ 


Note. — Average number of NR steps taken per solution, rate 
of full zone (cell) updates, and average failure rate (per solution) 
are shown for each method. The Intel C Compiler Version 8 and 
an Intel Xeon 3.06GHz workstation were used for these runs. 











Table 5. Comparison of Accretion Rates 


Method 

< Mo > 

< E/Mo > 

< L/Mo > 

2D 

-1.249 

0.86 

1.41 

IBw 

-1.145 

0.86 

1.32 


-1.235 

0.86 

1.43 

5D 

-1.226 

0.86 

1.40 

Thin disk 

— 

0.82 

1.95 


Table 6. Accretion Disk Efficiency Comparison 


Method 

NR steps per sol. 

Zone-cycles/node/sec. ^ 

Failure Rate 

2D 

4.19 

24535 

9.57 X 10"^ 

IBw 

4.18 

23860 

9.33 X 10"^ 


5.22 

20585 

9.46 X 10"^ 

5D 

4.52 

14741 

9.22 X 10"^ 


^Four nodes in parallel were used per run. The rates assume that run¬ 
time scales linearly with the number of nodes. 


Note. — Average number of NR steps taken per solution, rate of full 
zone (cell) per processor updates, and average failure rate (per solution) 
are shown for each method. The Intel C Compiler Version 7.1 and four 
Intel Xeon 2.40GHz processors were used for these runs. 
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10-10 10-8 10-6 10-4 10-2 lOO 


u 

Fig. 1.— Shown are the relative errors in calculating u from our parameter space survey averaged 
over {po, 7 , 4>} and plotted versus logio('*^)- Only those points for which all methods found a 
solution were included in the average, accounting for approximately 58% of the surveyed points. The 
curves corresponding to the 2D, IDw, 10*2, 10^2, 5D, and polynomial methods are represented 
by, respectively, circles, (blue in the electronic edition) exes, (red) squares, (magenta) asterisks, 
(green) empty triangles, and (cyan) filled triangles; please note that the 2D, ID^, 10*2, and 5D 
relative errors are almost indistinguishable. See the electronic edition of the Journal for a color 
version of this figure. 
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Fig. 2.— Histograms of the final value of FInr per parameter space point for the methods using a 
NR scheme. The distributions generated from the 2D, IDw, 10^2, lD.y 2 and 5D methods are rep¬ 
resented, respectively, by a solid line, (blue in the electronic edition) dots, (red) dashes, (magenta) 
long dashes, and (green) dot-dashes. Only those parameter space points for which all methods 
converged to a solution were included in this figure, so all points he below FInr = 10“^^. There are 
approximately 3.3 x 10® points per histogram. The first bin contains those points that lie beneath 
the displayed range. See the electronic edition of the Journal for a color version of this figure. 
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Fig. 3.— Histograms of the number of NR iterations taken by the methods per parameter space 
point. The distributions generated from the 2D, IDiy, 1D*2, 1D„2 and 5D methods are represented, 
respectively, by a solid line, (blue in the electronic edition) dots, (red) dashes, (magenta) long 
dashes, and (green) dot-dashes. Only those parameter space points for which all methods converged 
to a solution were included in this figure. There are approximately 3.3 x 10® points per histogram. 
See the electronic edition of the Journal for a color version of this figure. 




























- 29 - 




Fig. 4.— Snapshots of (top left), By (top right), logj^Qp (bottom left) and 7 (bottom right) 
taken at t = 4 from the cylindrical explosion evolution. A continuous greyscale from white to 
black is used for each plot, using the same limits as used by Komissarov (1999): Bx G [0.008,0.35], 
By G [—0.18,0.18], log^Qp G [—4.5, —1.5], 7 G [1,4.57]. The lower-right plot also shows the magnetic 
field lines that originate from x = —6 and along the y-axis at equal intervals of dy = 1/3 from 
-5.67 to 5.67. 
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Fig. 5.— Accretion rates of rest-mass, specific energy and specific angular momentum over time 
for a disk evolution around a Kerr black hole of spin parameter a = 0.9375M at a resolution of 
256^. The solid curve, (blue in the electronic version) dots, (red) dashes, and (green) dot-dashes 
represent runs that used—respectively—the 2D, IDw, 10^2 and 5D methods. All the straight 
horizontal lines indicate the time-averages of the accretion rates over t = 500 — 2000 GMjc? for the 
different methods except for the (cyan) lines of dots and long dashes which denote the thin disk 
values. See the electronic edition of the Journal for a color version of this figure. 






















Fig. 6.— Snapshots of log po sXt = 800 GM/c^ (top row) and t = 2000 GM/c^ (bottom row) on a 
256^ grid are shown from accretion disk evolutions using different methods: (from left to right) the 
2D, IDwi lh )*2 and 5D methods. The data are plotted here in {r,6) coordinates in the xy-plane. 
The region we exclude about the origin to excise the singularity from the grid can be seen near the 
origin. The color map is logarithmically-spaced such that the black (dark blue in the electronic 
version) points refer to po — 4 x 10“'^, while the white (dark red in the electronic version) regions 
correspond to po — 0.69. See the electronic edition of the Journal for a color version of this figure. 
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Fig. 7.— Histograms of the values of F'nr with which the methods return, over the entire spacetime 
region covered by the accretion disk evolution at a resolution of 128^. The solid curve, (blue 
in the electronic version) dots, (red) dashes, and (green) dot-dashes represent runs that used— 
respectively—the 2D, ID^y, 1D*2 and 5D methods. Again, points with F^nr > are considered 

erroneous and are replaced by interpolated values during the simulation (see text for further details). 
Approximately 3.8 x 10® points are represented here. The leftmost bin contains those points that 
lie below the displayed range. 
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Fig. 8.— Histograms of the number of NR iterations taken by the methods over the entire spacetime 
region covered in the accretion disk evolution at a resolution of 128^. Approximately 3.8 x 10® points 
are represented here. The bin located at Nnr = 31 contains all those points for which no solution 
was found, all other points were successful. 




































